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Abstract 

The multilayer multiconhguration time-dependent Hartree theory within second quantization 
representation of the Fock space is applied to study correlated electron transport in models of single- 
molecule junctions. Extending previous work, we consider models which include both electron- 
electron and electronic-vibrational interaction. The results show the influence of the interactions 
on the transient and the stationary electrical current. The underlying physical mechanisms are 
analyzed in conjunction with the nonequilibrium electronic population of the molecular bridge. 



I. INTRODUCTION 



The process of charge transport in molecular junctions has received much attention 
recently.— - — Single molecule junctions, consisting of single molecules that are chemically 
bound to metal electrodes, are well-suited systems to study nonequilibrium transport phe- 
nomena at the nanoscale and are also of interest for potential applications in the field of 
molecular electronics. Recent developments in experimental techniques, such as electromi- 
gration, mechanically controllable break junctions, or scanning tunneling microscopy,- 1 ^ - — 
have made it possible to study transport properties of molecular junctions. The rich 
experimental observations, e.g., Coulomb blockade,— Kondo effect,— negative differential 
resistance,— switching and hysteresis,— - — have stimulated many theoretical develop- 
ments for understanding quantum transport at the molecular scale. 

A particular challenge for the theory of charge transport in molecular junctions is the 
accurate treatment of correlation effects beyond the mean-field level. In molecular junctions, 
there are two types of correlation effects due to electronic-vibrational and electron-electron 
interaction. Considering vibrational induced correlation effects, a variety of theoretical ap- 
proaches have been developed, including scattering theory— - — nonequilibrium Green's func- 
tion approaches,— - — and master equation methods.— ^~— In spite of the physical insight 
offered by these methods, all of them involve significant approximations. For example, 
NEGF methods and master equation approaches are usually based on (self-consistent) per- 
turbation theory and/or employ factorization schemes. Scattering theory approaches to 
vibrationally coupled electron transport, on the other hand, neglect vibrational nonequilib- 
rium effects and are limited to the treatment of a small number of vibrational degrees of 
freedom. These shortcomings have motivated us to develop a systematic, numerically exact 
methodology to study quantum dynamics and quantum transport including many-body ef- 
fects — the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) theory in 
second quantization representation (SQR).— For a generic model of vibrationally coupled 
electron transport, we have demonstrated the importance of treating electronic-vibrational 
coupling accurately. Comparison with approximate methods such as NEGF reveals the ne- 
cessity of employing accurate methods such as the ML-MCTDH-SQR, in particular in the 
strong coupling regime. 

In this paper, we extend the ML-MCTDH-SQR method to treat electron-electron inter- 



action. Considering the paradigmatic Anderson impurity model, we show the applicability 
of the methodology to obtain an accurate description. Furthermore, we consider a model 
which incorporates both electron-electron and electronic- vibrational interaction. To the best 
of our knowledge, the results reported for this model are the first obtained by a numerically 
exact method. 

It is noted that a variety of other powerful methods have been developed in the recent 
years with the same goal, i.e., to facilitate numerically exact simulations for nonequilib- 
rium transport in model systems. These include the numerical path integral approach,—"— 
real-time quantum Monte Carlo simulations,— 1 ^ the numerical renormalization group 
approach,—, the time-dependent density matrix renormalization group approach.—, and 
the hierarchical equations of motion method^ 1 ^. For a comparison and an comprehen- 
sive overview of various different methods in the case of nonequilibrium transport with 



electron-electron interaction see Ref. 



75 



The remaining part of the paper is organized as follows. Section [TT] outlines the physical 
model and the observables of interest. The ML-MCTDH-SQR theory is described in Sec- 
tion IIII1 Section IIVI presents numerical results for a variety of parameter regimes as well as 
an analysis of the transport mechanisms. Section |V] concludes with a summary. 



II. MODEL AND OBSERVABLES OF INTEREST 



To study correlated electron transport in molecular junctions, we consider a generic model 
which includes both electron-electron and electronic-vibrational interaction. The model 
comprises two discrete electronic states (spin up and down) at the molecular bridge, two 
electronic continua describing the left and the right metal leads, respectively, and a distri- 
bution of harmonic oscillators that models the vibrational modes of the molecular bridge. 
The Hamiltonian reads 

H — H e i + H nvic + iJ e i_ nuc , (2.1a) 



where H e \, H nvLC , and H e \-nuc describe the electronic degrees of freedom, the nuclear vibra- 
tions, and their coupling terms, respectively 

H e \ = E d h d ,g + U d hd^h d)i + ^ E kL hk L)(T + ^ E kR h kRy(J (2.1b) 

a k L ,cr k R ,a 

+ Vdk L (dtc kL ,a + ct L> Ja) + V dkR (d+C kR>a + C+ Rt J a ), 

k L ,a k R ,a 
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H c i- nuc = Y ^d,a Y 2c iQj- (2.ld) 

In the above expression n denotes the number operator, subscript u d" refers to the bridge 
state, u ki/kji" the states of the left/right metal leads, and "a =t?4" the two spin states. 
Operators d + /d, c\Jc kL ^ c\Jc kR are the fermionic creation/annihilation operators for the 
electronic states on the molecular bridge, the left and the right leads, respectively. The sec- 
ond term in ( I2.1bl) describes the on-site Coulomb repulsion of the electrons on the molecular 
bridge with electron-electron coupling strength U d . The energies of the electronic states in 
the leads, E kL , E kR , as well as the molecule-lead coupling parameters V dkL , V dkR are assumed 
to be independent of the spin polarization and are defined through the energy-dependent 
level width functions 

T L (E) = 2tt J2\y<ik L \ 2 S(E-E kL ) } T R (E) = 2n J2\ydk R \ 2 S(E-E kR ). (2.2) 

k L k R 

Without electronic-vibrational interaction the model introduced above reduces to the 
well known Anderson impurity model,— which has been investigated in great detail both in 
equilibrium and nonequilibrium.— Neglecting, on the other hand, electron-electron inter- 
action, it corresponds to the standard model of vibrationally coupled electron transport in 
molecular junctions, which has also been studied in great detail, mostly based on approx- 
imate methods. Recently a numerically exact treatment of the latter model (i.e. without 
electron-electron interaction) became possible using path integral techniques,— 1 ^ as well as 
the ML-MCTDH approach.— iZ&>Z2, To the best of our knowledge, the full model including 
electron-electron and electronic-vibrational interaction has so far not been considered with 
numerically exact methods. 

In principle, the parameters of the model can be obtained for a specific molecular junction 
employing first-principles electronic structure calculations.— In this paper, which focuses on 



the methodology and general transport properties, however, we will use a generic parame- 
terization. Employing a tight-binding model, the function T(E) is given as 



, ^JAff-E 2 \E\ < 2I/3J 

T(E) = i fiV He 1 1 " '^ e| , (2.3a) 

\E\ >2\p e \ 

r L (E) = r(E - fi L ), r R (E) = t(e - p R ), (2.3b) 

where (3 e and a e are nearest-neighbor couplings between two lead sites and between the lead 
and the bridge state, respectively. I.e., the width functions for the left and the right leads 
are obtained by shifting T(E) relative to the chemical potentials of the corresponding leads. 
We consider a simple model of two identical leads, in which the chemical potentials are given 
by 

fi L/R = E f ± V/2, (2.4) 

where V is the source-drain bias voltage and Ef the Fermi energy of the leads. Since only 
the difference E& — Ef is physically relevant, we set Ef = 0. 

Similarly, the frequencies Uj and electronic-nuclear coupling constants Cj of the vibrational 
modes of the molecular junctions are modeled by a spectral density function^ 1 ^ 

J J 

In this paper, the spectral density is chosen in Ohmic form with an exponential cutoff 

JoM = |owe- w / w % (2.6) 

where a is the dimensionless Kondo parameter. 

Both the electronic and the vibrational continua can be discretized by choosing a density 
of states p e {E) and a density of frequencies p(u) such that^^r— 

Ek dEp e (E) = k, \ Vdk \ 2 = T ^ Ek l k=l,...,N e , (2.7a) 

2np e (E k ) 

du JP (u)=j, l = 2JoM j = l,...,N b . (2.7b) 

Uj 7T P {Uj) 

where N e is the number of electronic states (for a single spin/single lead) and Nj, is the 
number of bath modes in the simulation. In this work, we choose a constant p e {E), i.e., an 



equidistant discretization of the interval [—2(3 e ,2(3 e ], to discretize the electronic continuum. 
For the vibrational bath, p(u) is chosen as 

p( w ) = ^±l e -M_ (28) 

Within a given time scale the numbers of electronic states and bath modes are systematically 
increased to reach converged results for the quantum dynamics in the extended condensed 
phase system. In our calculations we employ 80-500 states for each electronic lead, implying 
40-250 electrons per lead, and a bath with 100-400 modes. 

The observable of interest for studying transport through molecular junctions is the 
current for a given source-drain bias voltage, given by (in this paper we use atomic units 
where h = e = 1) 

h(f) = JJ^ = ^[^ H[ HMe-^}, (2.9a) 

In(t) = d ^ = ^ tr {pe""i[H, N R ]e~""} . (2.9b) 
Here, Ni/n(t) denotes the time-dependent charge in each lead, defined as 

N ( {t) = ^tr[pe i&t N c e- idt ], ( = L,R } (2.10) 

and N( = J2 k( a n^ a is the occupation number operator for the electrons in each lead 
(( = L,R). For Hamiltonian ( 12. ip the explicit expression for the current operator is given 

as 

i ( =t[H,N ( ]=tY,V dk (dtc koa -ct c J (7 ), ( = L,R. (2.11) 

In the expressions above, p denotes the initial density matrix representing a grand- 
canonical ensemble for each lead and a certain preparation for the bridge state 



p = p° d exp [~(3(H - p L N L - p R N R )j , (2.12a) 



H0 = Y1 E ^k L ,a + E ^k R ,, + ^ n ° uc - (2.12b) 

Here p° d is the initial reduced density matrix for the bridge state, which is usually chosen as 
a pure state representing an occupied or an empty bridge state, and H® uc defines the initial 
bath equilibrium distribution. 



Various initial states can be considered. For example, one may choose an initially unoc- 
cupied bridge state and the nuclear degrees of freedom equilibrated with this state, i.e. an 
unshifted bath of oscillators with H® nc as given in Eq. (12. left . On the other hand, one may 
also start with a fully occupied bridge state and a bath of oscillators in equilibrium with the 
occupied bridge state 

(2.13) 



Other initial states may also be prepared. The initial state may affect the transient dynam- 
ics profoundly. The dependence of the steady-state current on the initial density matrix is a 
more complex issue. Recent investigations for a model without electron-electron interaction 
seem to indicate that different initial states may lead to different (quasi) steady states,— iSSiSl 
although this has been debated.— Even without coupling to a vibrational bath, the initial 
bridge state population may still affect the final stationary current in a time-dependent 
simulation.—*^ For all results reported in this paper, our calculations show that the sta- 
tionary state is independent on the initial condition within the error bar of the simulation 
(which is estimated to be less than 10% relative error). Since different sets of initial con- 
ditions also affect the time scale at which the current I(t) reaches its stationary value, we 
typically choose initial conditions that are close to the final steady state, e.g., an unoccupied 
initial bridge state if its energy is higher than the Fermi level of the leads and an occupied 
bridge state otherwise. 

The transient behavior of the thus defined currents and Ii{t) is usually different. 

However, the long-time limits of Ifi(t) and /l(£), which define the stationary current, are 
the same. It is found that the average current 

I(t)= 1 -[I R (t) + I L (t)}, (2.14) 

provides better numerical convergence properties by minimizing the transient characteristic, 
and thus will be used in most calculations. 

In our simulations the continuous set of electronic states of the leads is represented by a 
finite number of states. The number of states required to properly describe the continuum 
limit depends on the time t. The situation is thus similar to that of a quantum reactive 
scattering calculation in the presence of a scattering continuum, where, with a finite number 
of basis functions, an appropriate absorbing boundary condition is added to mimic the 



correct outgoing Green's function.— — Employing the same strategy for the present problem, 
the regularized electric current is given by 

J reg = lim r dt^-e-*. (2.15) 
J dt 

The regularization parameter t] is similar (though not identical) to the formal convergence 
parameter in the definition of the Green's function in terms of the time evolution operator 



G{E + ) = lim (-z) / dte i{E+ir >- H)t . (2.16) 
v->o+ J 

In numerical calculations, r] is chosen in a similar way as the absorbing potential used in 
quantum scattering calculations.— ~— In particular, the parameter rj has to be large enough to 
accelerate the convergence but still sufficiently small in order not to affect the correct result. 
While in the reactive scattering calculation t] is often chosen to be coordinate dependent, in 
our simulation rj is chosen to be time dependent 

[ (t < t) 

V(t) = { 1 ' (2.17) 

[r) -(t-r)/t (t>r). 

Here r/ is a damping constant, r is a cutoff time beyond which a steady state charge flow 
is approximately reached. As the number of electronic states increases, one may choose a 
weaker damping strength rjo and/or longer cutoff time r. The former approaches zero and 
the latter approaches infinity for an infinite number of states. In practice, for the systems 
considered in this work, convergence can be reached with a reasonable number of electronic 
states in the range of 80-500, with a typical r = 30-80 fs (a smaller r for less number of 
states) and l/r] = 3-10 fs. 

To gain insight into the transport mechanisms, it is also useful to consider the population 
of the electronic states localized on the the molecular bridge, which is given by 

P d (t) = ^tr jpe^ £ ^ e "^} • ( 2 - 18 ) 

III. THE MULTILAYER MULTICONFIGURATION TIME-DEPENDENT 
HARTREE THEORY IN SECOND QUANTIZATION REPRESENTATION 



The time-dependent study of transport properties in the model introduced above requires 
a method that is able to describe many-body quantum dynamics in an accurate and efficient 



way. For this purpose we employ the recently proposed Multilayer Multiconfiguration Time- 
Dependent Hartree Theory in Second Quantization Representation (ML-MCTDH-SQR),— 
which allows a numerically exact treatment of the many-body problem. Here we give a brief 
outline of the method. 

A. Overview of the ML-MCTDH theory 

The ML-MCTDH theory^ is a rigorous variational method to propagate wave packets 
in complex systems with many degrees of freedom. In this approach the wave function is 
represented by a recursive, layered expansion, 

l*W> = EE " E A nn- 3 M fl \<Pj2®)> ( 3 - la ) 

Q(k) 

i*S?(*)> = EE-E B Zt, Q jt) n i^^)' ( 3 - lb ) 

M{n,q) 

i4r ) (0)=EE- e c^^jt) n ( 3 - ic ) 

ai Q2 a M( K ,q) 7=1 

where A^...^^), A (t), Caia 2 q ...a M(Ktq) {t) and so on are the expansion coefficients for 
the first, second, third, layers, respectively; \ip^{t)), \v\ K,q \t)) , |6^ 9 ' 7 (t)), are the 
"single particle" functions (SPFs) for the first, second, third, layers. In Eq. (13.1 aft , p 
denotes the number of single particle (SP) groups/subspaces for the first layer. Similarly, 
Q(k) in Eq. (I3.1bl) is the number of SP groups for the second layer that belongs to the 
ftth SP group in the first layer, i.e., there are a total of YT K =i Q( K ) second layer SP groups. 
Continuing along the multilayer hierarchy, M(k, q) in Eq. (I3.1cl) is the number of SP groups 
for the third layer that belongs to the qih SP group of the second layer and the nth SP group 
of the first layer, resulting in a total of J2l =1 Hqif M{k, q) third layer SP groups. Naturally, 
the size of the system that the ML-MCTDH theory can treat increases with the number of 
layers in the expansion. In principle, such a recursive expansion can be carried out to an 
arbitrary number of layers. The multilayer hierarchy is terminated at a particular level by 
expanding the SPFs in the deepest layer in terms of time-independent configurations, each 
of which may contain several Cartesian degrees of freedom. 



The variational parameters within the ML-MCTDH theoretical framework are dynami- 
cally optimized through the use of the Dirac-Frenkel variational principle 9 ^ 

(8*(t)\i^-H\*(t))=0, (3.2) 

which results in a set of coupled, nonlinear differential equations for the expansion coefficients 
for all layers.— For a iV-layer version of the ML-MCTDH theory there are N + 1 levels 
of expansion coefficients. In this sense the conventional wave packet propagation method is 
a "zero-layer" MCTDH approach. 

The introduction of this recursive, dynamically optimized layering scheme in the ML- 
MCTDH wavefunction provides more flexibility in the variational functional, which results 
in a tremendous gain in our ability to study large many-body quantum systems. During 
the past few years, significant progress has been made in further development of the theory 
to simulate quantum dynamics and nonlinear spectroscopy of ultrafast electron transfer 
reactions in condensed phases.—^— The theory has also been generalized to study heat 
transport in molecular junctions^ and to calculate rate constants for model proton transfer 
reactions in molecules in solution . 109 ' 110 Recent work of Manthe has introduced an even 
more adaptive formulation based on a layered correlation discrete variable representation 
(CDVR).iii^ 

B. Treating identical particles using the second quantization representation of 
Fock space 

To extend the original ML-MCTDH approach to systems of identical quantum particles 
requires a method that incorporates the exchange symmetry explicitly. This is because an 
ordinary Hartree product within the first quantized picture is only suitable to describe a 
configuration for a system of distinguishable particles. One strategy is to employ a properly 
symmetrized wave function in the first quantized framework, i.e., permanents in a bosonic 
case or Slater determinants in a fermionic case. This led to the MCTDHF approach^ - — for 
treating identical fermions and the MCTDHB approach^ for treating identical bosons as 
well as combinations thereof.— However, this wave function-based symmetrization is only 
applicable to the single layer MCTDH theory but is incompatible with the ML-MCTDH 
theory with more layers — there is no obvious analog of a multilayer Hartree configuration 



if permanents/determinants are used to represent the wave function. As a result, the ability 
to treat much larger quantum systems numerically exactly was severely limited. 

To overcome this limitation we proposed a novel approach^ that follows a fundamentally 
different route to tackle many-body quantum dynamics of indistinguishable particles — an 
operator-based method that employs the second quantization formalism of many-particle 
quantum theory. This differs from many previous methods where the second quantization 
formalism is only used as a convenient tool to derive intermediate expressions for the first 
quantized form. In the new approach the variation is carried out entirely in the abstract 
Fock space using the occupation number representation. Therefore, the burden of handling 
symmetries of identical particles in a numerical variational calculation is shifted completely 
from wave functions to the algebraic properties of operators. 

The major difference between the ML-MCTDH-SQR theory for identical fermions and 
the previous ML-MCTDH theory for distinguishable particles is the way how operators act. 
For example, in the second quantized form the fermionic creation/annihilation operators 
fulfill the anti-commutation relations 

{d P ,a^} = a P a^ + d^dp = 5p Q , {dj,, d^} = {d P , d Q } = 0. (3.3) 

The symmetry of identical particles is thus realized by enforcing such algebraic properties 
of the operators. This can be accomplished by introducing a permutation sign operator 
associated with each fermionic creation/annihilation operator, which incorporates the sign 
changes of the remaining spin orbitals in all the SPFs whose subspaces are prior to it.— 
For example, if a purely electronic problem is considered and only one layer is present, the 
overall wave function and the SPFs have the form 

L 

i*w> = EE--E^i*..^(oni^ ) (*)>' ( 3 - 4a ) 

2 m " 11 i 

i^: } (*)> = E B THm { S) = E E - E mm-kj, (3.4b) 

J" K =1 ni=0n 2 =0 n mK =0 

where rij = 0, 1 are the occupation numbers. A fermionic creation operator is actually 
implemented as 

(a^) + = (ll S^j (4 K) ) + , (3.5) 



where (^t = 1, 2, k — 1) is the permutation sign operator that accounts for permuting 
{a^) + from the first subspace all the way through to the «th subspace, and (a^) + is the 
reduced creation operator that only takes care of the fermionic anti-commutation relation 
in the nth subspace. The operator-based anti-commutation constraint (I3.3P results in the 
following operations 



& K) )M H M) = EE - E ^ 

ni=0ri2=0 n mK =0 



v-1 



n<-i 

.9=1 



= EE- E 



rti=0 rt2=0 



n< -i) 



.9=1 



Bn[ j n 2 ...n m St) M K) • | 1„> ■ K ) , 

(3.6a) 

££t...n m a) |n!)|n 2 )..> mM ). (3.6b) 



I.e., (a[f^) + not only creates a particle in the z/th spin orbital if it is vacant, but also affects the 
sign of each term in this SPF according to where u is located and what the occupations are 
prior to it. Furthermore, the permutation sign operators S^, \i — 1, 2, k — 1, incorporate 
the sign changes of the remaining spin orbitals in all the SPFs whose subspaces are prior to 
that of (ai K) ) + . 

Thus, the occupation number states in the ML-MCTDH-SQR theory are treated in the 
same way as the degrees of freedom in the original ML-MCTDH theory, except that the 
orderings of all the SP groups in all layers need to be recorded and maintained in later 
manipulations. More importantly, the equations of motion have the same form as in the 
original ML-MCTDH theory. The only difference is that for identical fermions each cre- 
ation/annihilation operator of the Hamiltonian is effectively a product of operators: a re- 
duced creation/annihilation operator that only acts on the bottom-layer SPFs for the Fock 
subspace it belongs to, and a series of permutation sign operators that accounts for the 
fermionic anti-commutation relations of all the spin orbitals prior to it. In the multilayer 
case the implementation is sophisticated but can still be reduced to handling (many and 
complicated) basic building blocks in the MCTDH or ML-MCTDH theory — products of 
operators. Thereby, the action of each Hamiltonian term (product of creation/annihilation 
operators) can be split into a series of operations on individual Fock subspaces.— On the 
other hand, for identical bosons the implementation is much simpler because there is no sign 
change upon permutation. 

In the second quantized form, the wave function is represented in the abstract Fock 
space employing the occupation number basis. As a result, it can be expanded in the same 



multilayer form as that for systems of distinguishable particles. It is thus possible to extend 
the numerically exact treatment to much larger systems. The symmetry of the wave function 
in the first quantized form is shifted to the operator algebra in the second quantized form. 
The key point is that, for both phenomeno logical models and more fundamental theories, 
there are only a limited number of combination of fundamental operators. For example, in 
electronic structure theory only one- and two-electron operators are present. This means 
that one never needs to handle all, redundant possibilities of operator combinations as offered 
by the determinant form in the first quantized framework. It is exactly this property that 
provides the flexibility of representing the wave functions in multilayer form and treat them 
accurately and efficiently within the ML-MCTDH-SQR theory. It is also noted that the 
ML-MCTDH-SQR approach outlined above for fermions has also be formulated for bosons 
or combinations of fermions, bosons and distinguishable particles.— 

IV. RESULTS AND DISCUSSION 

In this section, we present applications of the ML-MCTDH-SQR methodology to the 
study of correlated electron transport employing the model described in Sec. [Ill In particular, 
we discuss the influence of electron-electron and electronic-vibrational interaction on the 
transport characteristics for selected examples. Unlike the noninteracting transport model 
(Ud = 0,Cj = 0), these results represent nontrivial solutions to a many-body quantum 
problem and are often beyond the perturbation treatment. All calculations presented in this 
paper are for zero temperature, which corresponds to the deep quantum regime, and is often 
the most challenging physical regime of the problem. Meanwhile, this regime is relatively 
easy for our approach since only one initial wave function is required. An investigation of 
systems at finite temperature as well as an analysis of the physical mechanisms in a broader 
parameter range will be the topic of future work. 

A. Effect of electron-electron interaction on transport characteristics for fixed 
nuclei 

We first focus on the influence of electron-electron interaction and consider models with- 
out electronic-vibrational coupling (cj = 0), i.e. for fixed nuclei. Fig. [TJ shows the time- 



dependent current and the corresponding bridge state population for a model with the 
following set of electronic parameters: The tight-binding parameters for the function T(E) 
are a e = 0.2 eV, (3 e = 1 eV, corresponding to a moderate molecule-lead coupling and a 
bandwidth of 4 eV. The energy of the bridge states is located 0.5 eV above Fermi energy 
and the source-drain voltage is V — 0.1 V, i.e. the model is in the off-resonant transport 
regime. The results for both the current and the population show pronounced transient 
oscillations that decay on a time scale of ~ T _1 and approach a stationary plateau at longer 
times, which represents the steady state. The overall values of the current and population are 
rather small because the transport takes place in the off-resonant regime. The comparison of 
the results obtained for different parameters Ud shows that for this model electron-electron 
interaction has no significant influence on the population and the current, and this includes 
both the transient behavior and the long-time stationary value. Qualitatively this can be 
understood from the fact that the model is in the off-resonant transport regime. At zero 
coupling strength Ud-, the bare energies of the electronic bridge states are the same and are 
outside the conductance window defined by the chemical potentials of the two electrodes. 
Including the on-site repulsion term removes the degeneracy of these two bridge states if 
one state is occupied. That is, when one of the bridge states is populated, the electronic 
energy of the other state is increased by the value of Ud- However, due to the fact that the 
initial bridge states are relatively far away from the conductance window, their populations 
are small. As a result, the overall electronic correlation effect is small for this set of model 
parameters. At a finer scale it can be seen that with the increase in Ud both the stationary 
current and bridge state population decreases. This is consistent with the fact that upon 
increasing Ud the energy of the doubly occupied state {Ed + Ud) is moved to higher energies 
and thus even further away from the conduction window. 

Figure |2] shows the time-dependent current and the corresponding bridge state population 
for another model, where the parameters are the same as in Fig. [1] except for Ed — Ef = 0, 
i.e., energy of the bridge states is located at the Fermi energy of the leads. For zero on-site 
coupling strength (Ud = 0), this set of parameters corresponds to the resonant tunneling 
regime and involves mixed electron/hole transport. This results in a significantly larger 
stationary current and a population of approximately one, because each bridge state has 
50% probability to be occupied. In this parameter regime, electron-electron interaction has 
a pronounced influence on the transport characteristics. Upon increase of Ud the steady 



state value of both the current and the population decreases significantly. This is due to 
the fact that for increasing Ud the energy of the doubly occupied state moves out of the 
conductance window. For interaction strengths Ud 3> O.leV, the bridge state can only be 
singly occupied, resulting in an overall population of n d = 1/2, and the doubly occupied 
state does not contribute to the current. 

We next consider in Figs. 131 IH a model with the same parameters as in the previous two 
cases except that the energy of the bridge state is below the Fermi energy, E d — Ef = —0.5 
eV. For vanishing electron-electron interaction, Ud = 0, this is again a non-resonant case. 
However, because the bridge state is, in contrast to Fig. [TJ located below the Fermi energy it 
is almost doubly occupied when Ud = 0. While the stationary current for Ud = (full black 
line in Fig. [3(a) ) is, due to particle- hole symmetry, in fact identical to that in Fig. Ufa), 
the dependence of the transport characteristics on the electron-electron interaction is more 
complex than in the above two cases. Upon moderate increase of Ud, the energy of the 
doubly occupied state, Ed + Ud, moves closer to the Fermi energy and enters the conductance 
window. As a result, the population of this state decreases as shown in Fig. [3](b). The 
current (Fig. [3(a)), on the other hand, increases because the doubly occupied state provides 
a channel for resonant transport. It is also observed that upon moderate increase of Ud 
the transient dynamics undergoes a coherent to incoherent transition. When Ud is further 
increased (Ud > 0.6 eV), the energy of the doubly occupied state becomes higher than 
the chemical potentials of both electrodes. As shown in Fig. HJ this causes a decrease of 
the current and the population. For large values of Ud, the population in the steady state 
approaches a value of unity, because the bridge state can only be singly occupied. It is 
interesting to note that for large Coulomb repulsion, e.g. Ud = 2eV, the initial transient 
current is negative. This is because in the simulation the bridge state is initially fully 
occupied. During the early transient time, electrons flow from the bridge states to both the 
left and the right electrodes, resulting in a negative transient current. As the bridge states 
approach their steady state population, electrons move continuously from the left electrode 
to the right electrode, establishing a steady-state current. 



B. Aspects of Coulomb blockade 



An interesting many-body nonequilibrium effect in charge transport in mesoscopic and 
nanosystems is Coulomb blockade.— ~— This phenomenon involves the suppression of the 
electrical current due to electron-electron interaction. Within the single-site Anderson im- 
purity model, the underlying mechanism is that the Coulomb repulsion with an electron 
that already occupies the bridge state prevents a second electron to transfer onto the bridge 
and thus reduces the current compared to a noninteracting model. 

This basic aspect of Coulomb blockade is demonstrated in Fig. [5j which shows simulated 
current- voltage characteristics for a resonant transport model, where the energy of the bridge 
states Ed is at the Fermi energy of the leads, Ed — Ef = 0. The tight-binding parameters 
for the function T(E) are a e — 0.1 eV, (3 e = 1 eV, corresponding to a smaller molecule- lead 
coupling and a bandwidth of 4 eV. Besides the noninteracting model (Ud = 0), three values 
of the electron-electron coupling strengths are considered: Ud = 0.5, 1, and 4 eV. To obtain 
the current- volt age characteristics, the stationary plateau value from the time- dependent 
simulation of the current was taken for each given voltage. The results show that upon 
inclusion of electron-electron interaction, the currents are suppressed at all voltages. The 
ratio between the blocked and unblocked currents attain a stationary ratio of approximately 
2/3 in the plateau region (within the convergence range of less than 10% relative error), and 
is nearly independent of the electronic coupling strength Ud- Within a zeroth order picture, 
this result can be rationalized as follows. 

For the model without electron-electron interaction {Ud = 0), there are three channels 
for electron transport through the two bridge states: (i) Electron transport through an 
initially unoccupied state. There are two such channels corresponding to the two spin 
polarizations of the bridge state, (ii) Electron transport through an initially singly occupied 
bridge resulting in the third channel which involves double occupation of the bridge state. 
When the source-drain voltage V is small, roughly \eV\ < 2Ud in the zeroth-order picture, 
the third two-electron transport channel is essentially closed, resulting in a current value 
2/3 of that for the unblocked case. At approximately \eV\ = 2Ud, e.g., V — 1 V for the 
case Ud = 0.5 eV in Fig. [51 the two-electron transport channel becomes available and the 
current begins to increase with the source-drain voltage. For finite molecule-lead coupling, 
the transition is broadened as shown in Fig. [51 For larger values of Ud, the energy of the 



doubly occupied state is outside the conductance window of the bias voltages considered 
and thus the current is suppressed. 

While Fig. [5] demonstrates the phenomenon of Coulomb blockade for varying the source- 
drain voltage, it is also instructive to study the phenomenon for varying the gate voltage. 
Assuming that an additional gate voltage V g predominantly shifts the energy of the bridge 
states E d , we can investigate the influence of the gate voltage by varying E d relative to 
the Fermi energy of the leads. The result depicted in Fig. exhibits the well known peak 
structures of Coulomb blockade, with maxima at energies E d = and E d = —U d , where 
the singly occupied levels and the doubly occupied level are in the conductance window, 
respectively. The parameters here are the same as in Fig. [5] except for a fixed U d = 0.5eV 
and a source-drain voltage V = 0.1V. It is noted that the value of the voltage considered is 
already beyond the linear response regime. 

C. Effect of electron-electron interaction on transport characteristics in the pres- 
ence of electron-vibrational coupling 

We finally consider a model which includes both electron-electron and electron-vibrational 
interaction. The presence of both interactions increases the complexity significantly. To the 
best of our knowledge, the results presented here are the first numerical exact simulations 
for this type of models. 

Fig. [7J shows results for a model, where the electronic parameters are the same as in 
Fig. El i.e., a e = 0.2 eV, (3 e = 1 eV, E d — Ef = 0, and a source-drain voltage of V = 0.1 V. 
The electronic degrees of freedom are coupled to a vibrational bath modeled by an Ohmic 
spectral density, as described Sec. [XTJ The characteristic frequency and the reorganization 
energy of the vibrational bath are u c = 500 cm -1 and A = 2au c = 0.25eV, respectively. 
These values are typical for larger molecular systems. Without Coulomb repulsion and 
coupling to the vibrational bath, this model corresponds to the resonant transport regime. 
Including the couplings to the vibrational modes has a significant impact on the electrical 
current. After a short transient time the coupling to the vibrations becomes effective and 
results in a suppression of the current. As illustrated by the solid black line in Fig. [7J the 
effect is very pronounced and the stationary current is essentially blocked. The underlying 
mechanism can be qualitatively rationalized by considering the energy level of the bridge 



states. For any finite bias voltage, the bare energy of the bridge states (E d — Ef = 0) is 
located between the chemical potential of the leads and thus, within a purely electronic 
model, current can flow. The coupling to the vibrations results in a polaron shift of the 
energy of the bridge state given by the reorganization energy A. For electronic- vibrational 
coupling strengths with A > \V\/2 the polaron-shifted energy of the bridge state is below 
the chemical potentials of both leads and thus current is blocked. This effect, referred to as 
phonon blockade of the current, has been observed e.g. in quantum dots^ and has been 
analyzed previously— As shown in Fig. EJb), the bridge states are almost fully occupied in 
this case. 

When the Coulomb repulsion term is included in the simulation (in addition to the 
vibrational bath), the energy level of the doubly occupied bridge state is shifted to higher 
energies as discussed for the previous models. For smaller values of Ud, this brings the 
polaron-shifted bridge state back to the conduction window and thus increases the stationary 
current. This can be seen from the currents for U d = 0.5eV and U d = leV in Fig. [7(a) and 
in Fig. [BJ which shows the current-voltage characteristics for a few selected values of Ud- It 
is evident that for small Ud the stationary current increases versus Ud- 

However, if Ud becomes too large, e.g., Ud = 2eV in Fig. Eta), the doubly occupied 
bridge state has too high energy and is located above the conduction window, which again 
results in a suppression of the stationary current. On the other hand, the overall popula- 
tion of the bridge state decreases monotonically upon increase of Ud and reaches a value 
of unity for large because then the bridge state can only be singly occupied. Due to 
the strongly correlated dynamics in this parameter regime, including both electron-electron 
and electronic-vibrational coupling, convergence for larger values of Ud and larger voltages 
than those depicted in Figs, [7J El is difficult within our present implementation of the 
ML-MCTDH-SQR methodology. Experience shows that convergence in this regime can 
be facilitated by transforming the current Hamiltonian to another form in order to reduce 
correlation effects. This will be the subject of future work. 

Although the interpretation of the above vibronic and electronic correlated transport 
properties is appealing in terms of the energetics of the bridge states, it should be empha- 
sized that the mechanism involves the formation of correlated many-body states that are 
significantly more complex than this noninteracting electronic picture, and cannot be fully 
described by just considering the static shift of the energy of the bridge states. This is 



evident by examining the strength of the interaction parameters A and Ud in Fig. [7J Thus, 
an accurate description of the vibrational and electronic dynamics as well as their couplings 
is essential to obtain a quantitative description of the many-body quantum dynamics and 
the transport characteristics. 

V. CONCLUDING REMARKS 

In this paper we have employed the ML-MCTDH-SQR method to study correlated elec- 
tron transport through model single-molecule junctions. Extending our previous work,— ^ 
we have considered models which include both electron-electron and electron-vibrational 
interaction. The ML-MCTDH-SQR method allows an accurate, in principle numerically 
exact treatment of this many-body quantum transport problem including both the transient 
dynamics and the steady state. 

The results obtained for selected model systems demonstrate the complex interplay of 
electronic and vibrational dynamics. For example, strong electron-vibrational coupling may 
result in a pronounced suppression of the electrical current (phonon blockade), which is 
accompanied by the formation of a polaron-like state. Including electron-electron interaction, 
this suppression of the current can be partially lifted because the transport channel provided 
by the doubly occupied bridge state shifts into the conductance window. 

In the present work we have considered a model with a single electronic state at the 
molecular bridge. It should be noted, however, that the ML-MCTDH-SQR method can also 
be applied to more complex models with various electronic states and interacting electrons in 
the leads. In addition to transport problems it may also be used to describe photoinduced 
dynamics in molecular adsorbates at metal or semiconductor surfaces including a proper 
description of correlation effects. Another important phenomenon in correlated electron 
transport is the Kondo effect.— iZLi22 The application of the methodology to simulate trans- 
port in the Kondo regime, in particular for very small voltage, requires special discretization 
techniques (e.g., the scheme pioneered by Wilson 123 ) and can be facilitated by the use of 
correlated initial states. This will be considered in future work. 
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FIG. 1: (a) Time-dependent current I(t) for different electron-electron coupling strength Ud and 
(b) the corresponding electronic population at the bridge state. Other parameters are: a e = 0.2eV, 
j3 e = leV, Ed — Ef = 0.5eV, and the source-drain voltage V = 0.1V. The inset in panel (a) depicts 
the stationary current in an enlarged view. 





FIG. 2: Same as Fig. Q] except for E d - E f = 0. 



(a) 



(b) 




40 

Time (fs) 




— u d 


= 


-- u d 


= 0.1 eV 


— u d 


= 0.2 eV 


... u d 


= 0.3 eV 


— u d 


= 0.4 eV 


— u d 


= 0.5 eV 


— u d 


= 0.6 eV 



40 

Time (fs) 



60 



FIG. 3: Same as Fig. Q] except for E d -E f = -0.5eV. 
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FIG. 4: Same as Fig. [3] except for a larger range of Ud- 
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FIG. 5: Current-voltage characteristics for different electron-electron coupling strength Ud- Other 
parameters are: a e = O.leV, (3 e = leV, E^ — Ej = 0. The lines are intended as a guide to the eye. 




FIG. 6: Dependence of the current on the gate voltage. The electronic parameters are: a e = O.leV, 
(3 e = leV, Ud = 0.5eV, and — Ef = for zero gate voltage. The source-drain voltage is 0.1V. 
The line is intended as a guide to the eye. 
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FIG. 7: (a) Time-dependent current I(t) for different electron-electron coupling strength and 
(b) the corresponding electronic population of the bridge state. The results are obtained for a 
model, which includes both electron-electron and electron-vibrational coupling. For comparison, 
result for a purely electronic model (i.e. Ud = 0, A = 0) are shown as indicated in the legend. 
The source-drain voltage is V = 0.1V. The electronic parameters are are: a e = 0.2eV, /3 e = leV, 
Ed — Ef = 0. The reorganization energy and characteristic frequency for the vibrational bath are 
A = 0.25 eV and u c = 500cm _1 , respectively. 
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FIG. 8: Current-voltage characteristics for different electron-electron coupling strength Ud- The 
results are obtained for a model, which includes both electron-electron and electron-vibrational 
coupling. Other parameters are: a e = 0.2eV, /3 e = leV, Ed — Et = 0. The reorganization 



energy and characteristic frequency for the vibrational bath are A = 0.25 eV and lo c 
respectively. The lines are intended as a guide to the eye. 
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